The database used here is composed of 39 Assaf ewes in the first lactation which had the estrus synchronized and the lambing performed in a short window in order to avoid variations caused by the lactation stage. These animals were housed in individual pens and were milked twice a day in a 1x10 stall milking parlor (DeLaval). In the first stage, these animals were fed ad libitum a total mixed ratio (TMR) formulated from alfalfa hay (particle size > 4 cm) and concentrates (50:50 forage:concentrate ratio). These ewes had their dry matter intake (DMI) and milk yield measured daily for a period of three weeks after a period of three weeks of adaptation to the TMR. For each animal, the feed intake was estimated daily by weighting the dry matter refused. Additionally, morning and evening milk production was weighted for each animal in order to calculate the daily milk yield. Protein, fat, and lactose content were estimated for each animal as described by Barrio et al. (2022). The changes in the body weight (BW) were calculated by each ewe through the recording of two consecutive days at the beginning and at the end of the FE evaluation.
#Replace the pathway by the folder in your computer
FE.db<-read.table("~/post_doc_Leon/projetos/RFI_sheep_R_calculations/databases_FE/FE_AllMetrics_database_sheep_07_02_23.txt", h=T, sep="\t")
head(FE.db)
## Animal_ID Treatment meanBW BWC1 days_period Fat Protein Milk intake
## 1 61472 Control 66.9 -1.0 31 66.2 49.1 2.63 2.610
## 2 61473 Control 73.4 0.0 31 59.8 46.1 2.55 2.714
## 3 61474 Control 66.1 0.6 31 55.3 49.8 1.53 2.462
## 4 61477 Challenge 54.4 -6.0 31 71.7 50.6 2.18 2.307
## 5 61480 Control 66.8 0.0 31 55.4 49.6 2.63 2.985
## 6 61483 Control 77.9 0.0 31 52.8 49.6 3.04 3.156
## energyTMR DIM
## 1 0.93 23
## 2 0.93 29
## 3 0.93 32
## 4 0.93 31
## 5 0.93 50
## 6 0.93 28
The following columns are observed in the database:
Animal_ID: Animal identification
Treatment: Categorical variable indicating if the animals was subjected or not to a protein restriction (~40% reduction) challenge
meanBW: Average body weight (kg)
BWC1: Body weight change during the experimental period
days_period: Days in the experiment
Fat: Fat yeild (g/kg)
Protein: Protein yield (g/kg)
Milk: Milk yield (kg/d)
intake: Actual dry matter intake
energyTMR: Energy content of the total mixed ratio (TMR)
DIM: Days in milk
library(ggplot2)
library(PerformanceAnalytics)
library(plotly)
library(patchwork)
#Defining custom theme for plots in the document
custom_theme<-function(){
theme(panel.grid.major = element_line(colour = "#d3d3d3"),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 14, face = "bold"),
axis.title = element_text(face="bold", size=20),
axis.text.x = element_text(colour="black", size = 20, angle = 0),
axis.text.y = element_text(colour="black", size = 20),
axis.line = element_line(size=0.5, colour = "black"),
legend.title = element_text(size = 25),
legend.text = element_text(size=20),
strip.text.x = element_text(size = 15))
}
In order to calculate the FE metrics, some additional values must be estimated.
Body weight change during the period (BWC2)
#Defined as the body weight change during the whole period divided by the days
#in the period
FE.db$BWC2<-FE.db$BWC1/FE.db$days_period
Metabolic body weight
#Defined as the body weight raised to 0.75
FE.db$MBW<-FE.db$meanBW^0.75
Energy requirements for BW change, UFL/d
#Defined based on INRA equations
FE.db$ER_BWC<-(0.28*(FE.db$BWC2*1000))/100
Energy requirements for maintenance, UFL/d
#Defined based on INRA equations
FE.db$ER_main<-(0.0345*FE.db$MBW)
Mean milk yield, L/d (MY)
FE.db$MY<-FE.db$Milk/1.036
Fat yield, L/d (MY)
FE.db$FY<-FE.db$Fat/1.036
Protein yield, L/d (MY)
FE.db$PY<-FE.db$Protein/1.036
Metabolic body weight change per day
#Defined as the absolute BWC1 raised to 0.75
FE.db$MBWC1<-abs(FE.db$BWC1)^0.75
Metabolic body weight change for the experimental period
#Defined as the absolute BWC2 raised to 0.75
FE.db$MBWC2<-abs(FE.db$BWC2)^0.75
Energy corrected milk (ECM, kg/d) - INRA 2018
#Defined based on INRA equations
FE.db$ECM<-FE.db$MY*(0.0071*FE.db$FY+0.0043*FE.db$PY+0.2224)
Energy requirements for milk yield, UFL/d
#Defined based on INRA equations
FE.db$ER_MY<-0.686*FE.db$ECM
Total energy requirements, UFL/d
#Defined as the sum of the Energy requirements for MY, maintenance, and BWC
FE.db$TER<-rowSums(FE.db[,c("ER_MY","ER_main","ER_BWC")])
Predicted DM intake
#Defined as the ratio between total energy requirements and the energy on TMR
FE.db$Pred_DMI<-FE.db$TER/FE.db$energyTMR
Checking new database
head(FE.db)
## Animal_ID Treatment meanBW BWC1 days_period Fat Protein Milk intake
## 1 61472 Control 66.9 -1.0 31 66.2 49.1 2.63 2.610
## 2 61473 Control 73.4 0.0 31 59.8 46.1 2.55 2.714
## 3 61474 Control 66.1 0.6 31 55.3 49.8 1.53 2.462
## 4 61477 Challenge 54.4 -6.0 31 71.7 50.6 2.18 2.307
## 5 61480 Control 66.8 0.0 31 55.4 49.6 2.63 2.985
## 6 61483 Control 77.9 0.0 31 52.8 49.6 3.04 3.156
## energyTMR DIM BWC2 MBW ER_BWC ER_main MY FY
## 1 0.93 23 -0.03225806 23.39212 -0.09032258 0.8070281 2.538610 63.89961
## 2 0.93 29 0.00000000 25.07680 0.00000000 0.8651495 2.461390 57.72201
## 3 0.93 32 0.01935484 23.18201 0.05419355 0.7997794 1.476834 53.37838
## 4 0.93 31 -0.19354839 20.03084 -0.54193548 0.6910640 2.104247 69.20849
## 5 0.93 50 0.00000000 23.36589 0.00000000 0.8061232 2.538610 53.47490
## 6 0.93 28 0.00000000 26.22123 0.00000000 0.9046325 2.934363 50.96525
## PY MBWC1 MBWC2 ECM ER_MY TER Pred_DMI
## 1 47.39382 1.0000000 0.07611649 2.233674 1.5323003 2.249006 2.418286
## 2 44.49807 0.0000000 0.00000000 2.027122 1.3906056 2.255755 2.425543
## 3 48.06950 0.6817316 0.05189102 1.193408 0.8186778 1.672651 1.798549
## 4 48.84170 3.8336586 0.29180462 1.943903 1.3335172 1.482646 1.594243
## 5 47.87645 0.0000000 0.00000000 2.051046 1.4070175 2.213141 2.379721
## 6 47.87645 0.0000000 0.00000000 2.318505 1.5904942 2.495127 2.682932
After the estimation of the additional metrics above mentioned, the estimation of FCR and FEI it is a straightforward approach.
#FEI is defined as Actual dry matter intake subtracted by the predicted dry matter intake
FE.db$FEI<-FE.db$intake-FE.db$Pred_DMI
#Checking distribution of FEI values
ggplot(FE.db, aes(x=FEI)) +
geom_histogram(color="black", fill="white") +
custom_theme()
#Checking correlation between FEI and actual intake values
cor_FEI_DMI<-round(cor(FE.db$intake, FE.db$FEI),2)
text_FEI<-paste("Animal ID: ", FE.db$Animal_ID, "\n","DMI: ",
round(FE.db$intake,2), "\n" ,
"FEI: ", round(FE.db$FEI,2),sep="")
FEI_plot<-ggplot(FE.db, aes(y=intake,x=FEI, text=text_FEI)) +
geom_point() +
custom_theme() + scale_y_continuous(name="DMI") +
ggtitle(paste("Actual DMI vs FEI (r2= ", cor_FEI_DMI,")",sep=""))
ggplotly(FEI_plot, tooltip = "text")
#FCR is defined as Feed intake (kg DM/d) / ECM (kg/d), where ECM is estimated based on the ECM equation from INRA 2018 (for sheep)
FE.db$FCR<-FE.db$intake/FE.db$ECM
#Checking distribution of FCR values
ggplot(FE.db, aes(x=FCR)) +
geom_histogram(color="black", fill="white") +
custom_theme()
#Checking correlation between FCR and actual intake values
cor_FCR_DMI<-round(cor(FE.db$intake, FE.db$FCR),2)
text_FCR<-paste("Animal ID: ", FE.db$Animal_ID, "\n","DMI: ",
round(FE.db$intake,2), "\n" ,
"FCR: ", round(FE.db$FCR,2),sep="")
FCR_plot<-ggplot(FE.db, aes(y=intake,x=FCR, text=text_FCR)) +
geom_point() +
custom_theme() + scale_y_continuous(name="DMI") +
ggtitle(paste("Actual DMI vs FCR (r2= ", cor_FCR_DMI,")",sep=""))
ggplotly(FCR_plot, tooltip = "text")
The RFI is defined as the subtraction of the actual feed intake by the predicted feed intake. The prediction of feed intake can be performed using different models. Here, for models to predict the feed intake and subsequently the RFI will be shown.
RFI estimation based on the model proposed by Pryce et al. 2015
#A linear model including ECM, MBW, BWC1 and DIM is used to predict the feed intake
model.pryce<-lm(intake ~ ECM + MBW + BWC1 + DIM,data=FE.db)
#Checking model output
summary(model.pryce)
##
## Call:
## lm(formula = intake ~ ECM + MBW + BWC1 + DIM, data = FE.db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35819 -0.08819 -0.00099 0.09000 0.26783
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.159225 0.347147 0.459 0.64939
## ECM 0.503667 0.057933 8.694 3.70e-10 ***
## MBW 0.061756 0.012793 4.827 2.87e-05 ***
## BWC1 0.057814 0.017037 3.393 0.00177 **
## DIM 0.004807 0.003190 1.507 0.14102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1377 on 34 degrees of freedom
## Multiple R-squared: 0.8102, Adjusted R-squared: 0.7879
## F-statistic: 36.29 on 4 and 34 DF, p-value: 7.939e-12
#Extracting residuals (RFI)
FE.db$RFI_pryce<-residuals(model.pryce)
head(FE.db$RFI_pryce)
## [1] -0.17160200 -0.15426355 0.08155537 0.12954168 0.10939969 0.07510774
#Checking distribution of RFI_pryce values
ggplot(FE.db, aes(x=RFI_pryce)) +
geom_histogram(color="black", fill="white") +
custom_theme()
#Checking the correlation between predicted and actual values
r2_pryce<-round(summary(model.pryce)$adj.r.squared,3)
RFI_pryce_table<-data.frame(DMI=FE.db$intake, predicted_DMI=predict(model.pryce))
text_pryce<-paste("Animal ID: ", FE.db$Animal_ID, "\n","DMI: ",
round(RFI_pryce_table$DMI,2), "\n" ,
"RFI_pryce: ", round(RFI_pryce_table$predicted_DMI,2),sep="")
plot_pryce<-ggplot(RFI_pryce_table, aes(y=DMI,x=predicted_DMI, text=text_pryce)) +
geom_point() +
custom_theme() +
ggtitle(paste("Actual DMI vs predicted DMI by pryce (r2= ", r2_pryce,")",sep=""))
ggplotly(plot_pryce, tooltip = "text")
RFI estimation based on ECM, MBW and MBWC2
#A linear model including ECM, MBW and MBWC2 is used to predict the feed intake
model.RFI2<-lm(intake ~ ECM + MBW + MBWC2,data=FE.db)
#Checking model output
summary(model.RFI2)
##
## Call:
## lm(formula = intake ~ ECM + MBW + MBWC2, data = FE.db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52406 -0.05510 0.00208 0.09454 0.31614
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.40214 0.34773 1.156 0.255
## ECM 0.43302 0.07291 5.939 9.29e-07 ***
## MBW 0.06804 0.01484 4.584 5.60e-05 ***
## MBWC2 -0.15832 0.42091 -0.376 0.709
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1751 on 35 degrees of freedom
## Multiple R-squared: 0.684, Adjusted R-squared: 0.657
## F-statistic: 25.26 on 3 and 35 DF, p-value: 7.073e-09
#Extracting residuals (RFI)
FE.db$RFI2<-residuals(model.RFI2)
head(FE.db$RFI2)
## [1] -0.33895494 -0.27219170 -0.02603707 -0.25362495 0.10486107 -0.03423561
#Checking distribution of RFI_pryce values
ggplot(FE.db, aes(x=RFI2)) +
geom_histogram(color="black", fill="white") +
custom_theme()
#Checking the correlation between predicted and actual values
r2_RFI2<-round(summary(model.RFI2)$adj.r.squared,3)
RFI_RFI2_table<-data.frame(DMI=FE.db$intake, predicted_DMI=predict(model.RFI2))
text_RFI2<-paste("Animal ID: ", FE.db$Animal_ID, "\n","DMI: ",
round(RFI_RFI2_table$DMI,2), "\n" ,
"RFI2: ", round(RFI_RFI2_table$predicted_DMI,2),sep="")
plot_RFI2<-ggplot(RFI_RFI2_table, aes(y=DMI,x=predicted_DMI, text=text_RFI2)) +
geom_point() +
custom_theme() +
ggtitle(paste("Actual DMI vs predicted DMI by RFI2 (r2= ", r2_RFI2,")",sep=""))
ggplotly(plot_RFI2, tooltip = "text")
RFI estimation based on ECM, MBW and an interaction between MeanBW and BWC2
#A linear model including ECM, MBW and an interaction term between meanBW and BWC2 is used to predict the feed intake
model.RFI3<-lm(intake ~ ECM + MBW + meanBW*BWC2,data=FE.db)
#Checking model output
summary(model.RFI3)
##
## Call:
## lm(formula = intake ~ ECM + MBW + meanBW * BWC2, data = FE.db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33820 -0.09806 0.01144 0.07646 0.27749
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.62071 3.91880 0.414 0.682
## ECM 0.52318 0.06219 8.413 1.02e-09 ***
## MBW -0.14470 0.69578 -0.208 0.837
## meanBW 0.05270 0.18517 0.285 0.778
## BWC2 -0.62360 3.81764 -0.163 0.871
## meanBW:BWC2 0.04517 0.06165 0.733 0.469
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1427 on 33 degrees of freedom
## Multiple R-squared: 0.8022, Adjusted R-squared: 0.7723
## F-statistic: 26.77 on 5 and 33 DF, p-value: 1.004e-10
#Extracting residuals (RFI)
FE.db$RFI3<-residuals(model.RFI3)
head(FE.db$RFI3)
## [1] -0.24251485 -0.20657496 0.04239889 0.05597538 0.15213866 0.01144312
#Checking distribution of RFI_pryce values
ggplot(FE.db, aes(x=RFI3)) +
geom_histogram(color="black", fill="white") +
custom_theme()
#Checking the correlation between predicted and actual values
r2_RFI3<-round(summary(model.RFI3)$adj.r.squared,3)
RFI_RFI3_table<-data.frame(DMI=FE.db$intake, predicted_DMI=predict(model.RFI3))
text_RFI3<-paste("Animal ID: ", FE.db$Animal_ID, "\n","DMI: ",
round(RFI_RFI3_table$DMI,2), "\n" ,
"RFI3: ", round(RFI_RFI3_table$predicted_DMI,2),sep="")
plot_RFI3<-ggplot(RFI_RFI3_table, aes(y=DMI,x=predicted_DMI, text=text_RFI3)) +
geom_point() +
custom_theme() +
ggtitle(paste("Actual DMI vs predicted DMI by pryce (r2= ", r2_RFI3,")",sep=""))
ggplotly(plot_RFI3, tooltip = "text")
The Pearson correlation across the FE variables estimated here will be calculated and plotted using the R package PerformanceAnalytics. It is possible to note that, in general, a strong correlation is observed between the variables.
chart.Correlation(FE.db[,c("FEI","FCR", "RFI_pryce", "RFI2", "RFI3")], histogram=TRUE, pch=19)
Let’s create a final databse containing only the Animal Id and the FE metrics estimated in the last steps.
#Filtering the database created in this tutorial based on the column names
FE.metrics<-FE.db[,c("Animal_ID","Treatment","FEI","FCR", "RFI_pryce", "RFI2", "RFI3")]
head(FE.metrics)
## Animal_ID Treatment FEI FCR RFI_pryce RFI2 RFI3
## 1 61472 Control 0.1917141 1.168478 -0.17160200 -0.33895494 -0.24251485
## 2 61473 Control 0.2884569 1.338844 -0.15426355 -0.27219170 -0.20657496
## 3 61474 Control 0.6634509 2.063000 0.08155537 -0.02603707 0.04239889
## 4 61477 Challenge 0.7127573 1.186788 0.12954168 -0.25362495 0.05597538
## 5 61480 Control 0.6052787 1.455355 0.10939969 0.10486107 0.15213866
## 6 61483 Control 0.4730680 1.361222 0.07510774 -0.03423561 0.01144312
This database could be save in your computer using the following code:
#Filtering the database created in this tutorial based on the column names
write.table(FE.metrics, file="Feed_EfficiencyMetrics_smarter.txt", col.names=T, row.names=F, sep="\t", quote=F)
This will create the tab-delimited file “Feed_EfficiencyMetrics_smarter.txt” in your work directory.
The animals which composed the database analyzed here were subject to a nutritional protein restriction (NPR). Initially, 40 lamb ewes from a single flock in the northwest region of Castilla y León (Spain) that were transported to the facilities of the IGM in León were fed a standard diet for replacement ewe lambs providing 16% crude protein until three months of age and were subsequently divided into two groups. The two experimental groups were composed of 20 NPR and 20 C animals. To evaluate the impact of feed restriction challenge due to a trade market problem and a shortage of concentrate inputs, the C ewes received the standard diet mentioned above for 64 d; during the same period, the NPR ewes received the same diet without soybean meal (44% reduction in protein intake). The 64-d NPR period in the prepuberal stage was coincident with the allometric growth of the mammary gland. Now, we will evaluate if the feed efficiency metrics estimated here are significantly different between control and challenge groups.
# Here we will compare the means of FEI between control and challenge groups using a t-test. But first, we need to define if there is variance equality between groups.
var.test(FEI ~ Treatment, data=FE.metrics)
##
## F test to compare two variances
##
## data: FEI by Treatment
## F = 1.5248, num df = 19, denom df = 18, p-value = 0.3758
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5918188 3.8816358
## sample estimates:
## ratio of variances
## 1.524777
#The result results of the variance equality test (p-value = 0.3758) suggest that the variances are equal between groups. Therefore, we are assuming this in the t-test.
t.test(FEI ~ Treatment, data=FE.metrics, var.equal = TRUE)
##
## Two Sample t-test
##
## data: FEI by Treatment
## t = 1.4407, df = 37, p-value = 0.1581
## alternative hypothesis: true difference in means between group Challenge and group Control is not equal to 0
## 95 percent confidence interval:
## -0.02965722 0.17561563
## sample estimates:
## mean in group Challenge mean in group Control
## 0.4925271 0.4195479
# Here we will compare the means of FCR between control and challenge groups using a t-test. But first, we need to define if there is variance equality between groups.
var.test(FCR ~ Treatment, data=FE.metrics)
##
## F test to compare two variances
##
## data: FCR by Treatment
## F = 1.0131, num df = 19, denom df = 18, p-value = 0.9814
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3932081 2.5789829
## sample estimates:
## ratio of variances
## 1.013071
#The result results of the variance equality test (p-value = 0.9814) suggest that the variances are equal between groups. Therefore, we are assuming this in the t-test.
t.test(FCR ~ Treatment, data=FE.metrics, var.equal = TRUE)
##
## Two Sample t-test
##
## data: FCR by Treatment
## t = -0.14205, df = 37, p-value = 0.8878
## alternative hypothesis: true difference in means between group Challenge and group Control is not equal to 0
## 95 percent confidence interval:
## -0.1705026 0.1481624
## sample estimates:
## mean in group Challenge mean in group Control
## 1.451066 1.462236
# Here we will compare the means of RFI_pryce between control and challenge groups using a t-test. But first, we need to define if there is variance equality between groups.
var.test(RFI_pryce ~ Treatment, data=FE.metrics)
##
## F test to compare two variances
##
## data: RFI_pryce by Treatment
## F = 1.1357, num df = 19, denom df = 18, p-value = 0.7907
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.4408033 2.8911515
## sample estimates:
## ratio of variances
## 1.135697
#The result results of the variance equality test (p-value = 0.7907) suggest that the variances are equal between groups. Therefore, we are assuming this in the t-test.
t.test(RFI_pryce ~ Treatment, data=FE.metrics, var.equal = TRUE)
##
## Two Sample t-test
##
## data: RFI_pryce by Treatment
## t = 1.3434, df = 37, p-value = 0.1873
## alternative hypothesis: true difference in means between group Challenge and group Control is not equal to 0
## 95 percent confidence interval:
## -0.02819448 0.13913472
## sample estimates:
## mean in group Challenge mean in group Control
## 0.02702390 -0.02844622
# Here we will compare the means of RFI2 between control and challenge groups using a t-test. But first, we need to define if there is variance equality between groups.
var.test(RFI2 ~ Treatment, data=FE.metrics)
##
## F test to compare two variances
##
## data: RFI2 by Treatment
## F = 0.49525, num df = 19, denom df = 18, p-value = 0.1377
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1922247 1.2607679
## sample estimates:
## ratio of variances
## 0.4952525
#The result results of the variance equality test (p-value = 0.1377) suggest that the variances are equal between groups. Therefore, we are assuming this in the t-test.
t.test(RFI2 ~ Treatment, data=FE.metrics, var.equal = TRUE)
##
## Two Sample t-test
##
## data: RFI2 by Treatment
## t = 0.45002, df = 37, p-value = 0.6553
## alternative hypothesis: true difference in means between group Challenge and group Control is not equal to 0
## 95 percent confidence interval:
## -0.08576579 0.13474031
## sample estimates:
## mean in group Challenge mean in group Control
## 0.01192969 -0.01255757
# Here we will compare the means of RFI3 between control and challenge groups using a t-test. But first, we need to define if there is variance equality between groups.
var.test(RFI3 ~ Treatment, data=FE.metrics)
##
## F test to compare two variances
##
## data: RFI3 by Treatment
## F = 0.87033, num df = 19, denom df = 18, p-value = 0.765
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3378039 2.2155967
## sample estimates:
## ratio of variances
## 0.8703264
#The result results of the variance equality test (p-value = 0.765) suggest that the variances are equal between groups. Therefore, we are assuming this in the t-test.
t.test(RFI3 ~ Treatment, data=FE.metrics, var.equal = TRUE)
##
## Two Sample t-test
##
## data: RFI3 by Treatment
## t = 1.6375, df = 37, p-value = 0.11
## alternative hypothesis: true difference in means between group Challenge and group Control is not equal to 0
## 95 percent confidence interval:
## -0.01620261 0.15271109
## sample estimates:
## mean in group Challenge mean in group Control
## 0.03325207 -0.03500218
We can also check the distribution of each value within control and challenge groups using boxplots.
#FEI
FEI_NPR_plot<-ggplot(FE.metrics, aes(y=FEI,x=Treatment, fill=Treatment)) +
geom_boxplot() + scale_x_discrete(name = "NPR Groups") +
theme_bw() + scale_fill_manual(values=c("chartreuse", "brown3")) +
ggtitle("FEI")
#FCR
FCR_NPR_plot<-ggplot(FE.metrics, aes(y=FCR,x=Treatment, fill=Treatment)) +
geom_boxplot() + scale_x_discrete(name = "NPR Groups") +
theme_bw() + scale_fill_manual(values=c("chartreuse", "brown3")) +
ggtitle("FCR")
#RFI_pryce
RFI_pryce_NPR_plot<-ggplot(FE.metrics, aes(y=RFI_pryce,x=Treatment, fill=Treatment)) +
geom_boxplot() + scale_x_discrete(name = "NPR Groups") +
theme_bw() + scale_fill_manual(values=c("chartreuse", "brown3")) +
ggtitle("RFI_pryce")
#RFI2
RFI2_NPR_plot<-ggplot(FE.metrics, aes(y=RFI2,x=Treatment, fill=Treatment)) +
geom_boxplot() + scale_x_discrete(name = "NPR Groups") +
theme_bw() + scale_fill_manual(values=c("chartreuse", "brown3")) +
ggtitle("RFI2")
#RFI3
RFI3_NPR_plot<-ggplot(FE.metrics, aes(y=RFI3,x=Treatment, fill=Treatment)) +
geom_boxplot() + scale_x_discrete(name = "NPR Groups") +
theme_bw() + scale_fill_manual(values=c("chartreuse", "brown3")) +
ggtitle("RFI3")
(FEI_NPR_plot | FCR_NPR_plot) / (RFI_pryce_NPR_plot | RFI2_NPR_plot | RFI3_NPR_plot)
If we check the summary of the models used here to estimate RFI, we will notice that the intercept is not significant for any model. Consequently, what would happen if we exclude the intercept from our models? But, even more important… Should we remove the intercept from our model in this case?
RFI estimation based on the model proposed by Pryce et al. 2015
#A linear model including ECM, MBW, BWC1 and DIM is used to predict the feed intake
model.pryce_NoInt<-lm(intake ~ 0 + ECM + MBW + BWC1 + DIM,data=FE.db)
#Checking model output
summary(model.pryce_NoInt)
##
## Call:
## lm(formula = intake ~ 0 + ECM + MBW + BWC1 + DIM, data = FE.db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35945 -0.09059 0.00206 0.09164 0.27497
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ECM 0.503955 0.057272 8.799 2.16e-10 ***
## MBW 0.066926 0.005981 11.190 4.13e-13 ***
## BWC1 0.055159 0.015842 3.482 0.00136 **
## DIM 0.005724 0.002456 2.331 0.02567 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1361 on 35 degrees of freedom
## Multiple R-squared: 0.9979, Adjusted R-squared: 0.9977
## F-statistic: 4184 on 4 and 35 DF, p-value: < 2.2e-16
#Extracting residuals (RFI)
FE.db$RFI_pryce_NoInt<-residuals(model.pryce_NoInt)
head(FE.db$RFI_pryce_NoInt)
## [1] -0.15772512 -0.15188646 0.09281117 0.14026517 0.10135052 0.07240085
#Checking distribution of RFI_pryce_NoInt values
ggplot(FE.db, aes(x=RFI_pryce_NoInt)) +
geom_histogram(color="black", fill="white") +
custom_theme()
#Checking the correlation between predicted and actual values
r2_pryce_NoInt<-round(summary(model.pryce)$adj.r.squared,3)
RFI_pryce_NoInt_table<-data.frame(DMI=FE.db$intake, predicted_DMI=predict(model.pryce))
text_pryce_NoInt<-paste("Animal ID: ", FE.db$Animal_ID, "\n","DMI: ",
round(RFI_pryce_NoInt_table$DMI,2), "\n" ,
"RFI_pryce_NoInt: ",
round(RFI_pryce_NoInt_table$predicted_DMI,2),sep="")
plot_pryce_NoInt<-ggplot(RFI_pryce_NoInt_table, aes(y=DMI,x=predicted_DMI, text=text_pryce_NoInt)) +
geom_point() +
custom_theme() +
ggtitle(paste("Actual DMI vs predicted DMI by pryce (r2= ", r2_pryce_NoInt,")",sep=""))
ggplotly(plot_pryce_NoInt, tooltip = "text")
RFI estimation based on ECM, MBW and MBWC2
#A linear model including ECM, MBW and MBWC2 is used to predict the feed intake
model.RFI2_NoInt<-lm(intake ~ 0 + ECM + MBW + MBWC2,data=FE.db)
#Checking model output
summary(model.RFI2_NoInt)
##
## Call:
## lm(formula = intake ~ 0 + ECM + MBW + MBWC2, data = FE.db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53295 -0.06949 0.01668 0.08753 0.33648
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ECM 0.454906 0.070744 6.430 1.86e-07 ***
## MBW 0.083409 0.006646 12.551 1.04e-14 ***
## MBWC2 -0.074708 0.416595 -0.179 0.859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1759 on 36 degrees of freedom
## Multiple R-squared: 0.9964, Adjusted R-squared: 0.9961
## F-statistic: 3335 on 3 and 36 DF, p-value: < 2.2e-16
#Extracting residuals (RFI)
FE.db$RFI2_NoInt<-residuals(model.RFI2_NoInt)
head(FE.db$RFI2_NoInt)
## [1] -0.35153773 -0.29977962 -0.01059957 -0.22624461 0.10304205 -0.08578759
#Checking distribution of RFI_pryce values
ggplot(FE.db, aes(x=RFI2_NoInt)) +
geom_histogram(color="black", fill="white") +
custom_theme()
#Checking the correlation between predicted and actual values
r2_RFI2_NoInt<-round(summary(model.RFI2_NoInt)$adj.r.squared,3)
RFI_RFI2_NoInt_table<-data.frame(DMI=FE.db$intake, predicted_DMI=predict(model.RFI2_NoInt))
text_RFI2_NoInt<-paste("Animal ID: ", FE.db$Animal_ID, "\n","DMI: ",
round(RFI_RFI2_NoInt_table$DMI,2), "\n" ,
"RFI2_NoInt: ",
round(RFI_RFI2_NoInt_table$predicted_DMI,2),sep="")
plot_RFI2_NoInt<-ggplot(RFI_RFI2_NoInt_table, aes(y=DMI,x=predicted_DMI, text=text_RFI2_NoInt)) +
geom_point() +
custom_theme() +
ggtitle(paste("Actual DMI vs predicted DMI by RFI2 (r2= ", r2_RFI2_NoInt,")",sep=""))
ggplotly(plot_RFI2_NoInt, tooltip = "text")
RFI estimation based on ECM, MBW and an interaction between MeanBW and BWC2
#A linear model including ECM, MBW and an interaction term between meanBW and BWC2 is used to predict the feed intake
model.RFI3_NoInt<-lm(intake ~ 0 + ECM + MBW + meanBW*BWC2,data=FE.db)
#Checking model output
summary(model.RFI3_NoInt)
##
## Call:
## lm(formula = intake ~ 0 + ECM + MBW + meanBW * BWC2, data = FE.db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34078 -0.08942 0.01195 0.06968 0.27717
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ECM 0.53103 0.05850 9.078 1.31e-10 ***
## MBW 0.14260 0.03842 3.711 0.000734 ***
## meanBW -0.02368 0.01313 -1.804 0.080091 .
## BWC2 -1.03377 3.64135 -0.284 0.778210
## meanBW:BWC2 0.05177 0.05882 0.880 0.384963
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1409 on 34 degrees of freedom
## Multiple R-squared: 0.9978, Adjusted R-squared: 0.9975
## F-statistic: 3123 on 5 and 34 DF, p-value: < 2.2e-16
#Extracting residuals (RFI)
FE.db$RFI3_NoInt<-residuals(model.RFI3_NoInt)
head(FE.db$RFI3_NoInt)
## [1] -0.24907915 -0.20006074 0.04175125 0.05170519 0.14589767 0.03058491
#Checking distribution of RFI_pryce values
ggplot(FE.db, aes(x=RFI3_NoInt)) +
geom_histogram(color="black", fill="white") +
custom_theme()
#Checking the correlation between predicted and actual values
r2_RFI3_NoInt<-round(summary(model.RFI3_NoInt)$adj.r.squared,3)
RFI_RFI3_NoInt_table<-data.frame(DMI=FE.db$intake, predicted_DMI=predict(model.RFI3_NoInt))
text_RFI3_NoInt<-paste("Animal ID: ", FE.db$Animal_ID, "\n","DMI: ",
round(RFI_RFI3_NoInt_table$DMI,2), "\n" ,
"RFI3_NoInt: ",
round(RFI_RFI3_NoInt_table$predicted_DMI,2),sep="")
plot_RFI3_NoInt<-ggplot(RFI_RFI3_NoInt_table, aes(y=DMI,x=predicted_DMI, text=text_RFI3_NoInt)) +
geom_point() +
custom_theme() +
ggtitle(paste("Actual DMI vs predicted DMI by RFI3 (r2= ", r2_RFI3_NoInt,")",sep=""))
ggplotly(plot_RFI3_NoInt, tooltip = "text")
It is posisble to note that much higher adjusted R-squares were obtained for all the models. However, should we trust on those values?
Problems when you remove the intercept
Let’s use the models with higher R-adjusted with intercept and its version without intercept as example. First, we should check the distribution of the residuals obtained in both models.
# Checking the mean and standard deviation of the residuals obtained in each model
table.res<-data.frame(res_pryce=residuals(model.pryce), res_pryce_NoInt=residuals(model.pryce_NoInt))
head(table.res)
## res_pryce res_pryce_NoInt
## 1 -0.17160200 -0.15772512
## 2 -0.15426355 -0.15188646
## 3 0.08155537 0.09281117
## 4 0.12954168 0.14026517
## 5 0.10939969 0.10135052
## 6 0.07510774 0.07240085
| Mean | SD | |
|---|---|---|
| RFI Pryce | 0e+00 | 0.1302 |
| RFI Pryce without intercept | 6e-04 | 0.1306 |
res.pryce<-ggplot(model.pryce, aes(y=.resid,x=.fitted)) +
geom_point() +
custom_theme() +
ggtitle("RFI Pryce et al. (2015)") +
custom_theme() + scale_y_continuous(name="Residuals") +
scale_x_continuous(name="Fitted values")
res.pryce.NoInt<-ggplot(model.pryce_NoInt, aes(y=.resid,x=.fitted)) +
geom_point() +
custom_theme() +
ggtitle("RFI Pryce et al. (2015) without intercept") +
custom_theme() + scale_y_continuous(name="Residuals") +
scale_x_continuous(name="Fitted values")
res.pryce | res.pryce.NoInt
Conclusion about intercept removal
As we can note, the SD obtained for the residuals in both models are quite similar. However, the mean of the residuals in the model without intercept is not equal zero, which is a requirement for the kind of model used here. Additionally, when we remove the intercept, we force the program to set the this value to be equal zero. Consequently, assuming that our the intercepts of our regression lines pass through the zero, which is not the case of most of the models fitted with real data (in the case of a continuous variable). In summary, to run a model without intercept, in the current case, result in an inflation of the sums of squares (SS) for the model (SSmodel) and residuals (SSresidual), leading to the increase in the R-square values.